#include "header.h"

static int nG, nP[128], nT, dims;
static grid_ptr G, H, LFT, RGHT;	
static FILE *fhlog;  
static double dt, a[8][8], b[8];
static fftw_complex m0,m1,m2,m3,m4,m5,m6,m7,m8,m9;
static fftw_complex l0,l1,l2,l3,l4,l5,l6,l7,l8,l9;
static fftw_complex t0,t1,t2,t3,h0,h1,h2,p0,p1; 
//static HMDT xx,yy,zz,aa;
static double max_mr, hyp_v;
static double rr[7];
static double r0s, aux, tau_0;

void construct_grids(){
  int j;
  G = malloc(sizeof(type_grid));
  RGHT = G;
  fprintf(fhlog, "#Assembled subgrid %d\n", 0);
  G->right = NULL;
  for (j = 0; j < nG-1; j++) {
    G->left = malloc(sizeof(type_grid));
    H = G;
    G = G->left;
    G->right = H;  
    fprintf(fhlog, "#Assembled subgrid %d\n", j+1);
  }
  LFT = G;
  G->left = NULL;
}

void init_f(input_ptr p){
  int j,k,cumN;
  double sm_a, sm_at, t_c = 0.5, lam = 4; 
  double rs, y, yp, ypp;

  lam   = p->a; t_c = p->sigma; 
  sm_a  = pow(50./3.,0.2)*pow(lam*t_c, 0.4);
  sm_at = -0.4*pow(50.*(lam*lam)/3.,0.2)*pow(t_c,-0.6);
  rs = sm_a*(1. - 0.5*pow(sm_a/(lam*lam), 1./3));
  y   = pow(sm_a, -1.5)*sqrt(lam*lam*(1. - pow(rs/sm_a, 2.)));
  yp  = -pow(lam, 2.)*rs/(y*pow(sm_a, 5.));
  ypp = -pow(lam, 2.)/pow(sm_a, 5)*(1./y + (pow(lam,2)*rs*rs)/(pow(sm_a, 5)*pow(y,3)) );

  /*xx = mgl_create_data();
  yy = mgl_create_data();
  zz = mgl_create_data();
  aa = mgl_create_data();*/
  nG = p->n_grids; nT = 0; 
  for (j = 0; j < nG; j++) {
    nP[j] = p->n_points[j];
    nT = nT + nP[j] + 8;
  }
  /*mgl_data_create(xx,nT,1,1);
  mgl_data_create(yy,nT,1,1);
  mgl_data_create(zz,nT,1,1);
  mgl_data_create(aa,nT,1,1);*/
  fhlog = fopen("run.log","a");
  construct_grids(); G = LFT;
  p->left = LFT; p->right = RGHT;
  dt = (p->cfl)*pow(p->h,2);
  G->a = 0; 
  G->h = p->h;
  G->n = nP[0];
  G->psi = malloc(sizeof(fftw_complex)*((G->n)+18));  
  if (p->zakh) {
    G->dens = malloc(sizeof(fftw_complex)*((G->n)+18));  
    G->vel = malloc(sizeof(fftw_complex)*((G->n)+18));  
  }
  G->wrk = malloc(sizeof(fftw_complex)*((G->n)+18));
  G->wrk2 = malloc(sizeof(fftw_complex)*((G->n)+18));
  G->rhs = malloc(sizeof(fftw_complex)*((G->n)+18));
  G->rk  = malloc(sizeof(fftw_complex)*((G->n)+18));
  for (int jj = 0; jj < 7; jj++) {
    G->k[jj] = malloc(sizeof(fftw_complex)*((G->n)+18));
    if (p->zakh) {
      G->kzd[jj] = malloc(sizeof(fftw_complex)*((G->n)+18));
      G->kzv[jj] = malloc(sizeof(fftw_complex)*((G->n)+18));
    }
  }
  G->r   = malloc(sizeof(double)*((G->n)+18));
  G->b = G->a + nP[0]*(G->h);
  fprintf(fhlog,"#Total number of points %d\n", nT);
  printf("a  %e\tdel r  %e\tdel xi  %e\n",sm_a, sm_a-rs, 1 - rs/sm_a  ); 
  for (k = 0; k < (G->n)+18; k++) {
    G->r[k] = G->a + (k-9)*(G->h);
    if (p->ic_type == 0) {
      G->psi[k] = (p->a)*exp(-pow((G->r[k])/(p->sigma), 2));
      if (p->zakh) {
        G->dens[k] = cabs(G->psi[k])*cabs(G->psi[k]);
        G->vel[k] = 0;
      }
    }
    if (p->ic_type == 1) {
      if ((G->r[k]) < rs) {
        G->psi[k] = pow(sm_a,-1.5)*sqrt(lam*lam*(1 - pow((G->r[k])/sm_a, 2.) )); 
        G->psi[k] = (G->psi[k])*cexp(0.25I*(sm_at/sm_a)*pow(G->r[k],2) );
      }
      else {
	G->psi[k] = (y + ((G->r[k])-rs)*yp)*cexp(-pow((1.*(G->r[k])-rs), 2) );
	G->psi[k] = (G->psi[k])*cexp(0.25I*(sm_at/sm_a)*pow(G->r[k],2) );
      }  
    }
    max_mr = cabs(G->psi[9]);
  }
  cumN = (G->n)+18;
  fprintf(fhlog,"#[%8.3f;%8.3f]\th=%.4f\t%5d pnts\n",G->a,G->b,G->h,G->n);
  j = 0; G = G->right;
  while (G != NULL) {
    G->a = (G->left)->b;
    G->h = 2*((G->left)->h);
    G->n = nP[j+1]; 
    G->psi = malloc(sizeof(fftw_complex)*((G->n)+18)); 
    if (p->zakh) {
      G->dens = malloc(sizeof(fftw_complex)*((G->n)+18));  
      G->vel = malloc(sizeof(fftw_complex)*((G->n)+18));  
    }
    G->wrk = malloc(sizeof(fftw_complex)*((G->n)+18));
    G->wrk2 = malloc(sizeof(fftw_complex)*((G->n)+18));
    G->rhs = malloc(sizeof(fftw_complex)*((G->n)+18));
    G->rk  = malloc(sizeof(fftw_complex)*((G->n)+18));
    for (int jj = 0; jj < 7; jj++) {
      G->k[jj] = malloc(sizeof(fftw_complex)*((G->n)+18));
      if (p->zakh) {
        G->kzd[jj] = malloc(sizeof(fftw_complex)*((G->n)+18));
        G->kzv[jj] = malloc(sizeof(fftw_complex)*((G->n)+18));
      }
    }
    G->r   = malloc(sizeof(double)*((G->n)+18));  
    for (k = 0; k < (G->n)+18; k++) {
      G->r[k] = G->a + (k-9)*(G->h);
      if (p->ic_type == 0) {
        G->psi[k] = (p->a)*exp(-pow((G->r[k])/(p->sigma), 2)); 	 
      }
      if (p->ic_type == 1) {
        if ((G->r[k]) < rs) {
          G->psi[k] = pow(sm_a,-1.5)*sqrt(lam*lam*(1 - pow((G->r[k])/sm_a,2.))); 
          G->psi[k] = (G->psi[k])*cexp(0.25I*(sm_at/sm_a)*pow(G->r[k],2) );
        }
        else {
	  //G->psi[k] = 0;
	  G->psi[k] = (y + ((G->r[k])-rs)*yp)*cexp(-pow(1.*((G->r[k])-rs), 2) );
          G->psi[k] = (G->psi[k])*cexp(0.25I*(sm_at/sm_a)*pow(G->r[k],2) );
        }
      } 
      if (p->zakh) {
        G->dens[k] = cabs(G->psi[k])*cabs(G->psi[k]);
        G->vel[k] = 0;
      }
    }
    cumN = cumN + (G->n)+18;
    G->b = G->a + nP[j+1]*(G->h);
    j++;
    fprintf(fhlog,"#[%8.3f;%8.3f]\th=%.4f\t%5d pnts\n",G->a,G->b,G->h,G->n);
    G = G->right;  
  }
  G = LFT; 
  if (p->ic_type == 1) setup_ic(p);   
  fclose(fhlog);
  dims = p->dims; r0s = 1./(R0*R0);
  p->tau = 0; p->m_p = cabs(LFT->psi[4]);
  hyp_v = p->hyp_v;
  /* l0 =    35./65536; --------------------------   10th order interpolation coefficients 
  l1 =  -405./65536;
  l2 =   567./16384;
  l3 = -2205./16384;
  l4 = 19845./32768; */

  l0 = -3.537570592015982e-07;   //  ---------------------  20th order interpolation coefficients
  l1 =  7.512135198339820e-06;
  l2 = -7.662377902306616e-05;
  l3 =  5.010016320738941e-04;
  l4 = -2.368371351622045e-03;
  l5 =  8.684028289280832e-03;
  l6 = -2.605208486784250e-02;
  l7 =  6.773542065639049e-02;
  l8 = -1.693385516409762e-01;
  l9 =  6.209080226835795e-01;

  t0 = -2.441406250000000e-003;   //  ----------------------  8th order interpolation
  t1 =  2.392578125000000e-002;
  t2 = -1.196289062500000e-001;
  t3 =  5.981445312500000e-001;

  h0 =  1.171875000000000e-002;   //   ----------------------  6th order
  h1 = -9.765625000000000e-002;
  h2 =  5.859375000000000e-001;

  p0 = -6.250000000000000e-002;   //   ----------------------  4th order
  p1 =  5.625000000000000e-001;


  for (int l = 0; l < 8; l++ ) for (int m = 0; m < 8; m++ ) a[l][m] = 0;

  a[2][1] =   4./7;

  a[3][1] = 115./112;
  a[3][2] =  -5./16;

  a[4][1] = 589./630;
  a[4][2] =   5./18;
  a[4][3] = -16./45;

  a[5][1] = 229./1200 -  29./6000 *sqrt(5);
  a[5][2] = 119./240  - 187./1200 *sqrt(5);
  a[5][3] = -14./75   +  34./375  *sqrt(5);
  a[5][4] =   0.      -   3./100  *sqrt(5);

  a[6][1] =  71./2400 - 587./12000*sqrt(5);
  a[6][2] = 187./480  - 391./2400 *sqrt(5);
  a[6][3] = -38./75   +  26./375  *sqrt(5);
  a[6][4] =  27./80   -   3./400  *sqrt(5);
  a[6][5] =   1./4    +   1./4    *sqrt(5); 

  a[7][1] = -49./480  +  43./160  *sqrt(5);
  a[7][2] =-425./96   +  51./32   *sqrt(5);
  a[7][3] =  52./15   -   4./5    *sqrt(5);
  a[7][4] = -27./16   +   3./16   *sqrt(5);
  a[7][5] =   5./4    -   3./4    *sqrt(5);
  a[7][6] =   5./2    -   1./2    *sqrt(5);

  b[1] = 1./12; b[2] = 0.; b[3] = 0.; b[4] = 0.;
  b[5] = 5./12; b[6] = 5./12; b[7] = 1./12;
}
//========================================================
void derivative_1() {
  // compute derivative in interior points
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < G->n; j++ ) {
      m4 =    (G->psi[j+5] - G->psi[j+13])/280;
      m3 = -4*(G->psi[j+6] - G->psi[j+12])/105;
      m2 =    (G->psi[j+7] - G->psi[j+11])/5;
      m1 = -4*(G->psi[j+8] - G->psi[j+10])/5;
      G->wrk[j+9] = (m1+m2+m3+m4)/(G->h);
    }
    G = G->right;
  }
}
void derivative_2() {
  // compute derivative in interior points
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < G->n; j++ ) {
      m4 =     -(G->psi[j+5] + G->psi[j+13])/560;
      m3 =    8*(G->psi[j+6] + G->psi[j+12])/315;
      m2 =     -(G->psi[j+7] + G->psi[j+11])/5;
      m1 =    8*(G->psi[j+8] + G->psi[j+10])/5;
      m0 = -205*(G->psi[j+9])/72;  
      m3 = m4 + m3;
      m2 = m3 + m2;
      m1 = m2 + m1;
      G->wrk[j+9] = (m0+m1)/((G->h)*(G->h));
    }
    G = G->right;
  }
}

void derivative_4() {
  // compute derivative in interior points
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < G->n; j++ ) {
      m4 =    7*(G->psi[j+5] + G->psi[j+13])/240;
      m3 =   -2*(G->psi[j+6] + G->psi[j+12])/5;
      m2 =  169*(G->psi[j+7] + G->psi[j+11])/60;
      m1 = -122*(G->psi[j+8] + G->psi[j+10])/15;
      m0 =   91*(G->psi[j+9])/8;  
      m3 = m4 + m3;
      m2 = m3 + m2;
      m1 = m2 + m1;
      G->wrk[j+9] = (m0+m1)/(pow(G->h,4));
    }
    G = G->right;
  }
}

void derivative_8() {
  // compute derivative in interior points
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < G->n; j++ ) {
      m4 =      (G->psi[j+5] + G->psi[j+13]);
      m3 =  -8.*(G->psi[j+6] + G->psi[j+12]);
      m2 =  28.*(G->psi[j+7] + G->psi[j+11]);
      m1 = -56.*(G->psi[j+8] + G->psi[j+10]);
      m0 =  70.*(G->psi[j+9]);  
      m3 = m4 + m3;
      m2 = m3 + m2;
      m1 = m2 + m1;
      G->wrk[j+9] = (m0+m1)/(pow(G->h,8));
    }
    G = G->right;
  }
}
void derivative_16() {
  // compute derivative in interior points
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < G->n; j++ ) {
      m8 =  1.000e+0*(G->psi[j+1] + G->psi[j+17]);
      m7 = -1.600e+1*(G->psi[j+2] + G->psi[j+16]);
      m6 =  1.200e+2*(G->psi[j+3] + G->psi[j+15]);
      m5 = -5.600e+2*(G->psi[j+4] + G->psi[j+14]);
      m4 =  1.820e+3*(G->psi[j+5] + G->psi[j+13]);
      m3 = -4.368e+3*(G->psi[j+6] + G->psi[j+12]);
      m2 =  8.008e+3*(G->psi[j+7] + G->psi[j+11]);
      m1 = -1.144e+4*(G->psi[j+8] + G->psi[j+10]);
      m0 =  1.287e+4*(G->psi[j+9]);  
      m7 = m8 + m7;
      m6 = m7 + m6;
      m5 = m6 + m5;
      m4 = m5 + m4;
      m3 = m4 + m3;
      m2 = m3 + m2;
      m1 = m2 + m1;
      G->wrk[j+9] = (m0+m1)/(pow(G->h,16));
    }
    G = G->right;
  }
}

void rhs() {
  derivative_2();
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < (G->n); j++) { 
      G->rhs[j+9] = 1.I*((G->wrk[j+9]) + pow(cabs(G->psi[j+9]),2)*(G->psi[j+9]));
    }
    if ((G->left)==NULL) G->rhs[9] = 1.I*(dims*(G->wrk[9]) 
                                   + pow(cabs(G->psi[9]),2)*(G->psi[9]));
    G = G->right;
  }
  derivative_1();
  G = LFT;
  while (G!=NULL) {
    if ((G->left)!=NULL) {
      for (int j = 0; j < (G->n); j++) { 
        G->rhs[j+9] += 1.I*(dims-1.)*(G->wrk[j+9])/(G->r[j+9]);
      }
    } else {
      for (int j = 1; j < (G->n); j++) { 
        G->rhs[j+9] += 1.I*(dims-1.)*(G->wrk[j+9])/(G->r[j+9]);
      }
    }
    G = G->right;
  }

  // Experimental Hyper-Viscosity
  derivative_16();
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < (G->n); j++) { 
      G->rhs[j+9] += -1.*(hyp_v*pow(G->h,16))*(G->wrk[j+9]);
    }
    G = G->right;
  }

  
}

void rhs_zakh() {
  //unfinished <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  HEEEEEREEEE 
  derivative_2();
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < (G->n); j++) { 
      G->rhs[j+4]   = 1.I*((G->wrk[j+4]) + pow(cabs(G->psi[j+4]),2)*(G->psi[j+4]));
      G->rhszd[j+4] = G->vel[j+4];
      G->rhszv[j+4] = 0; // 
    }
    if ((G->left)==NULL) G->rhs[4] = 1.I*(dims*(G->wrk[4]) 
                                   + pow(cabs(G->psi[4]),2)*(G->psi[4]));
    G = G->right;
  }
  derivative_1();
  G = LFT;
  while (G!=NULL) {
    if ((G->left)!=NULL) {
      for (int j = 0; j < (G->n); j++) { 
        G->rhs[j+4] += 1.I*(dims-1.)*(G->wrk[j+4])/(G->r[j+4]);
      }
    } else {
      for (int j = 1; j < (G->n); j++) { 
        G->rhs[j+4] += 1.I*(dims-1.)*(G->wrk[j+4])/(G->r[j+4]);
      }
    }
    G = G->right;
  }
}


void evolve_erk6(input_ptr p) {
  rhs();
  G = LFT; 

  tau_0    = p->tau;
  //rr[0]    = r0s*pow(cabs(G->psi[4]), 2); 
  find_ss(p);
  rr[0]    = pow(p->L,-2);
  aux      = a[2][1]*rr[0];
  p->tau   = tau_0 + aux*dt;

  while (G!=NULL) {
    memcpy(G->wrk2, G->psi, sizeof(fftw_complex)*((G->n)+18));
    memcpy(G->k[0], G->rhs, sizeof(fftw_complex)*((G->n)+18));
    for (int j = 0; j < (G->n); j++) { 
      G->rk[j+9] = a[2][1]*(G->k[0][j+9]);     
      G->psi[j+9] = G->wrk2[j+9] + dt*(G->rk[j+9]);
    }
    G = G->right; 
  }
  update_boundaries();
  for (int m = 3; m < 8; m++) {
    rhs();
    G = LFT;

    //rr[m-2]   = r0s*pow(cabs(G->psi[4]), 2); aux = 0;
    find_ss(p); 
    rr[m-2]   = pow(p->L,-2); aux = 0;
    
    for (int l = 1; l < m; l++) aux += a[m][l]*rr[l-1]; 
    p->tau   = tau_0 + aux*dt;

    while (G!=NULL) {
      memcpy(G->k[m-2], G->rhs, sizeof(fftw_complex)*((G->n)+18));
      for (int j = 0; j < (G->n); j++) { 
        G->rk[j+9] = 0;
        for (int l = 1; l < m; l++) G->rk[j+9] += a[m][l]*(G->k[l-1][j+9]);
        G->psi[j+9] = G->wrk2[j+9] + dt*(G->rk[j+9]);
      }
      G = G->right; 
    }
    update_boundaries();
  }
  rhs();
  G = LFT;

  //rr[6]   = r0s*pow(cabs(G->psi[4]), 2); aux = 0;
  find_ss(p);
  rr[6]   = pow(p->L,-2); aux = 0;
  for (int l = 1; l < 8; l++) aux += b[l]*rr[l-1]; 
  p->tau   = tau_0 + aux*dt;

  while (G!=NULL) {
    memcpy(G->k[6], G->rhs, sizeof(fftw_complex)*((G->n)+18));
    for (int j = 0; j < (G->n); j++) { 
      G->rk[j+9] = 0;
      for (int l = 1; l < 8; l++) G->rk[j+9] += b[l]*(G->k[l-1][j+9]);
      G->psi[j+9] = G->wrk2[j+9] + dt*(G->rk[j+9]);
    }
    G = G->right;
  }
  update_boundaries();
  find_ss(p);
}

void find_ss(input_ptr p) {
  G = LFT;
  p->m_p  = cabs(G->psi[9]);
  m4 =     -(cabs(LFT->psi[5]) + cabs(LFT->psi[13]))/560;
  m3 =    8*(cabs(LFT->psi[6]) + cabs(LFT->psi[12]))/315;
  m2 =     -(cabs(LFT->psi[7]) + cabs(LFT->psi[11]))/5;
  m1 =    8*(cabs(LFT->psi[8]) + cabs(LFT->psi[10]))/5;
  m0 = -205*(cabs(LFT->psi[9]))/72;  
  m3 = m4 + m3;
  m2 = m3 + m2;
  m1 = m2 + m1;
  p->p_rr = (m0+m1)/((LFT->h)*(LFT->h));
  p->v    = 1./sqrt(1. + 2.*(p->p_rr)/pow(p->m_p, 3));
  p->L    = (p->v)/(p->m_p);
}

void find_beta(input_ptr p) {
  double b0, b1, b2, b3, b4, b5, b6;
  b6  =   -5.71952248774065;
  b5  =   64.58276074952414;
  b4  = -286.42012849730742;
  b3  =  610.28639150972106;
  b2  = -581.01029967009026;
  b1  =   96.56198202862554;
  b0  =  131.01891277343202;

  p->beta   = b0 + (p->v)*(b1 + (p->v)*(b2 + (p->v)*(b3 + (p->v)*(b4 + (p->v)*(b5 + (p->v)*b6) ) ) ) );
}


double regrid(input_ptr p) {
 G = LFT;
 /*if (creal( (G->psi[4])*conj(G->rk[4])*2*dt) > (p->mr_tol)) add_grid(p); */
 if (cabs(G->psi[9]) > 2.*max_mr) {
   add_grid(p); G = LFT;
   max_mr = cabs(G->psi[9]);
   printf("%e\n", max_mr);
 }
 return dt;
}

void evolve_erk4(input_ptr p) {
  rhs();
  G = LFT;
  while (G!=NULL) {
    memcpy(G->wrk2, G->psi, sizeof(fftw_complex)*((G->n)+8));
    memcpy(G->rk  , G->rhs, sizeof(fftw_complex)*((G->n)+8));
    for (int j = 0; j < (G->n); j++) { 
      G->psi[j+4] = G->wrk2[j+4] + 0.5*dt*(G->rhs[j+4]);
    }
    G = G->right;
  }
  update_boundaries();
  rhs();
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < (G->n); j++) { 
      G->rk[j+4]  = G->rk[j+4] + 2*(G->rhs[j+4]);
      G->psi[j+4] = G->wrk2[j+4] + 0.5*dt*(G->rhs[j+4]);
    }
    G = G->right;
  }
  update_boundaries();
  rhs();
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < (G->n); j++) { 
      G->rk[j+4]  = G->rk[j+4] + 2*(G->rhs[j+4]);
      G->psi[j+4] = G->wrk2[j+4] + dt*(G->rhs[j+4]);
    }
    G = G->right;
  }
  update_boundaries();
  rhs();
  G = LFT;
  while (G!=NULL) {
    for (int j = 0; j < (G->n); j++) { 
      G->rk[j+4]  = G->rk[j+4] + (G->rhs[j+4]);
      G->psi[j+4] = G->wrk2[j+4] + (dt/6)*(G->rk[j+4]);
    }
    G = G->right;
  }
  update_boundaries();
}

void update_boundaries() {
  G = LFT;
  while (G!=NULL) {
    if (((G->left )==NULL) && ((G->right)==NULL)) {
      for (int j = 0; j < 9; j++) G->psi[j] = G->psi[18-j];
      G->psi[G->n+9] = G->psi[G->n+7];
      G->psi[G->n+10] = G->psi[G->n+6];
      G->psi[G->n+11] = G->psi[G->n+5];
      G->psi[G->n+12] = G->psi[G->n+4];
      G->psi[G->n+13] = G->psi[G->n+3];
      G->psi[G->n+14] = G->psi[G->n+2];
      G->psi[G->n+15] = G->psi[G->n+1];
      G->psi[G->n+16] = G->psi[G->n+0];
      G->psi[G->n+17] = G->psi[G->n-1];

    } else if (((G->left )==NULL) && ((G->right)!=NULL)) {
      for (int j = 0; j < 9; j++) G->psi[j] = G->psi[18-j];
      G->psi[G->n+9]  = (G->right)->psi[9];
      G->psi[G->n+11] = (G->right)->psi[10];
      G->psi[G->n+13] = (G->right)->psi[11];
      G->psi[G->n+15] = (G->right)->psi[12];
      G->psi[G->n+17] = (G->right)->psi[13];

      for ( int j = 0; j < 4; j ++) {
        /*m0 = l0*(G->psi[G->n-9+2*j] + (G->right)->psi[19+j]);  
        m1 = l1*(G->psi[G->n-7+2*j] + (G->right)->psi[18+j]);
        m2 = l2*(G->psi[G->n-5+2*j] + (G->right)->psi[17+j]);
        m3 = l3*(G->psi[G->n-3+2*j] + (G->right)->psi[16+j]);
        m4 = l4*(G->psi[G->n-1+2*j] + (G->right)->psi[15+j]);
        m5 = l5*(G->psi[G->n+1+2*j] + (G->right)->psi[14+j]);
        m6 = l6*(G->psi[G->n+3+2*j] + (G->right)->psi[13+j]);
        m7 = l7*(G->psi[G->n+5+2*j] + (G->right)->psi[12+j]);
        m8 = l8*(G->psi[G->n+7+2*j] + (G->right)->psi[11+j]);
        m9 = l9*(G->psi[G->n+9+2*j] + (G->right)->psi[10+j]);*/
        
        m0 = t0*(G->psi[G->n+3+2*j] + (G->right)->psi[13+j]);
        m1 = t1*(G->psi[G->n+5+2*j] + (G->right)->psi[12+j]);
        m2 = t2*(G->psi[G->n+7+2*j] + (G->right)->psi[11+j]);
        m3 = t3*(G->psi[G->n+9+2*j] + (G->right)->psi[10+j]);

        /*m0 = h0*(G->psi[G->n+3+2*j] + (G->right)->psi[13+j]);
        m1 = h1*(G->psi[G->n+5+2*j] + (G->right)->psi[12+j]);
        m2 = h2*(G->psi[G->n+7+2*j] + (G->right)->psi[11+j]);*/

        /*m0 = p0*(G->psi[G->n+3+2*j] + (G->right)->psi[13+j]);
        m1 = p1*(G->psi[G->n+5+2*j] + (G->right)->psi[12+j]);*/

        /*m1 = m1 + m0; m2 = m2 + m1; m3 = m3 + m2; m4 = m4 + m3; m5 = m5 + m4; 
        m6 = m6 + m5; m7 = m7 + m6; m8 = m8 + m7; */
        G->psi[G->n+10+2*j] = m0 + m1 + m2 + m3; 
      }

    } else if (((G->left )!=NULL) && ((G->right)==NULL)) {
      G->psi[G->n+9] = G->psi[G->n+7];
      G->psi[G->n+10] = G->psi[G->n+6];
      G->psi[G->n+11] = G->psi[G->n+5];
      G->psi[G->n+12] = G->psi[G->n+4];
      G->psi[G->n+13] = G->psi[G->n+3];
      G->psi[G->n+14] = G->psi[G->n+2];
      G->psi[G->n+15] = G->psi[G->n+1];
      G->psi[G->n+16] = G->psi[G->n+0];
      G->psi[G->n+17] = G->psi[G->n-1];

      G->psi[0] = (G->left)->psi[(G->left)->n - 9];
      G->psi[1] = (G->left)->psi[(G->left)->n - 7];
      G->psi[2] = (G->left)->psi[(G->left)->n - 5];
      G->psi[3] = (G->left)->psi[(G->left)->n - 3];
      G->psi[4] = (G->left)->psi[(G->left)->n - 1]; 
      G->psi[5] = (G->left)->psi[(G->left)->n + 1];
      G->psi[6] = (G->left)->psi[(G->left)->n + 3];
      G->psi[7] = (G->left)->psi[(G->left)->n + 5];
      G->psi[8] = (G->left)->psi[(G->left)->n + 7];     

    } else if (((G->left )!=NULL) && ((G->right)!=NULL)) {
      G->psi[G->n+9]  = (G->right)->psi[9];
      G->psi[G->n+11] = (G->right)->psi[10];
      G->psi[G->n+13] = (G->right)->psi[11];
      G->psi[G->n+15] = (G->right)->psi[12];
      G->psi[G->n+17] = (G->right)->psi[13];

      for ( int j = 0; j < 4; j ++) {
        /*m0 = l0*(G->psi[G->n-9+2*j] + (G->right)->psi[19+j]);  
        m1 = l1*(G->psi[G->n-7+2*j] + (G->right)->psi[18+j]);
        m2 = l2*(G->psi[G->n-5+2*j] + (G->right)->psi[17+j]);
        m3 = l3*(G->psi[G->n-3+2*j] + (G->right)->psi[16+j]);
        m4 = l4*(G->psi[G->n-1+2*j] + (G->right)->psi[15+j]);  20-th order (terrible)
        m5 = l5*(G->psi[G->n+1+2*j] + (G->right)->psi[14+j]);
        m6 = l6*(G->psi[G->n+3+2*j] + (G->right)->psi[13+j]);
        m7 = l7*(G->psi[G->n+5+2*j] + (G->right)->psi[12+j]);
        m8 = l8*(G->psi[G->n+7+2*j] + (G->right)->psi[11+j]);
        m9 = l9*(G->psi[G->n+9+2*j] + (G->right)->psi[10+j]);*/

        m0 = t0*(G->psi[G->n+3+2*j] + (G->right)->psi[13+j]);
        m1 = t1*(G->psi[G->n+5+2*j] + (G->right)->psi[12+j]);    //8-th order (best)
        m2 = t2*(G->psi[G->n+7+2*j] + (G->right)->psi[11+j]);
        m3 = t3*(G->psi[G->n+9+2*j] + (G->right)->psi[10+j]);

        /*m0 = h0*(G->psi[G->n+3+2*j] + (G->right)->psi[13+j]);
        m1 = h1*(G->psi[G->n+5+2*j] + (G->right)->psi[12+j]);    6-th order (worse)
        m2 = h2*(G->psi[G->n+7+2*j] + (G->right)->psi[11+j]);*/

        /*m0 = p0*(G->psi[G->n+3+2*j] + (G->right)->psi[13+j]);
        m1 = p1*(G->psi[G->n+5+2*j] + (G->right)->psi[12+j]);*/

        /*m1 = m1 + m0; m2 = m2 + m1; m3 = m3 + m2; m4 = m4 + m3; m5 = m5 + m4; 
        m6 = m6 + m5; m7 = m7 + m6; m8 = m8 + m7; */
        G->psi[G->n+10+2*j] = m0 + m1 + m2 + m3; 
      }


      /*
      m0 = l0*(G->psi[G->n-4] + (G->right)->psi[9]);
      m1 = l1*(G->psi[G->n-2] + (G->right)->psi[8]);
      m2 = l2*(G->psi[G->n+0] + (G->right)->psi[7]);
      m3 = l3*(G->psi[G->n+2] + (G->right)->psi[6]);
      m4 = l4*(G->psi[G->n+4] + (G->right)->psi[5]);
      G->psi[G->n+5] = m0 + m1 + m2 + m3 + m4;
      m0 = l0*(G->psi[G->n-2] + (G->right)->psi[10]);
      m1 = l1*(G->psi[G->n+0] + (G->right)->psi[9]);
      m2 = l2*(G->psi[G->n+2] + (G->right)->psi[8]);
      m3 = l3*(G->psi[G->n+4] + (G->right)->psi[7]);
      m4 = l4*(G->psi[G->n+6] + (G->right)->psi[6]);
      G->psi[G->n+7] = m0 + m1 + m2 + m3 + m4;
      */

      G->psi[0] = (G->left)->psi[(G->left)->n - 9];
      G->psi[1] = (G->left)->psi[(G->left)->n - 7];
      G->psi[2] = (G->left)->psi[(G->left)->n - 5];
      G->psi[3] = (G->left)->psi[(G->left)->n - 3];
      G->psi[4] = (G->left)->psi[(G->left)->n - 1]; 
      G->psi[5] = (G->left)->psi[(G->left)->n + 1];
      G->psi[6] = (G->left)->psi[(G->left)->n + 3];
      G->psi[7] = (G->left)->psi[(G->left)->n + 5];
      G->psi[8] = (G->left)->psi[(G->left)->n + 7];  
    }
    G = G->right;
  }
}
//==============================================================
/*void draw_tf(char *fname) {
 double p_max;
 int cumN = 0;
 HMGL gr = mgl_create_graph(800,600);
 G = LFT; p_max = cabs(G->psi[4]);
 while (G != NULL) {
   for (register long j = 0; j < (G->n)+8; j++){
     mgl_data_set_value(xx, (G->r[j])*p_max, j+cumN, 0, 0);
     mgl_data_set_value(yy, creal(G->psi[j])/p_max, j+cumN, 0, 0); // /pmax
     mgl_data_set_value(zz, cimag(G->psi[j])/p_max, j+cumN, 0, 0);
     mgl_data_set_value(aa, cabs(G->psi[j])/p_max, j+cumN, 0, 0);
   }
   cumN = cumN + (G->n) + 8;
   G = G->right;
 }

 //mgl_subplot(gr,2,1,0,"");
 mgl_set_ranges(gr, -0.5, 24.6, -1.2, 1.2, 0, 0);
 mgl_axis(gr,"xyz","","");
 mgl_box(gr);
 mgl_plot_xy(gr, xx, yy, "+b1", "");
 mgl_plot_xy(gr, xx, zz, "+g1", "");
 mgl_plot_xy(gr, xx, aa, "+r1", "");

 mgl_write_png_solid(gr,fname,"");
 mgl_delete_graph(gr);
}*/

void write_tf(char* fname) {
 G = LFT;
 FILE *fh = fopen(fname,"w");
 fprintf(fh, "#1.r 2.psi\n\n");
 while (G != NULL) {
   for (int j = 0; j < (G->n)+18; j++){
     fprintf(fh, "%.14e\t%.14e\t%.14e\t%.14e\n", G->r[j], cabs(G->psi[j]), creal(G->psi[j]), cimag(G->psi[j]));
   }
   G = G->right;
   fprintf(fh,"\n");
 }
 fclose(fh); 
}

void setup_ic(input_ptr p){
}

//===============================================================

void add_grid(input_ptr p){
  G = malloc(sizeof(type_grid));
  H = malloc(sizeof(type_grid));
  G->left = NULL;      G->right = H;
  H->left = G;         H->right = LFT->right; 
  if ((H->right)!=NULL) (H->right)->left = H;

  G->n = LFT->n;       H->n = (G->n)/2;
  G->h = 0.5*(LFT->h); H->h = (G->h)*2;
  G->a = 0;            G->b = (G->a) + (G->n)*(G->h);
  H->a = G->b;         H->b = (H->a) + (H->n)*(H->h);
 
  G->psi  = malloc(sizeof(fftw_complex)*((G->n)+18));  
  G->wrk  = malloc(sizeof(fftw_complex)*((G->n)+18));
  G->wrk2 = malloc(sizeof(fftw_complex)*((G->n)+18));
  G->rhs  = malloc(sizeof(fftw_complex)*((G->n)+18));
  G->rk   = malloc(sizeof(fftw_complex)*((G->n)+18));

  H->psi  = malloc(sizeof(fftw_complex)*((H->n)+18));  
  H->wrk  = malloc(sizeof(fftw_complex)*((H->n)+18));
  H->wrk2 = malloc(sizeof(fftw_complex)*((H->n)+18));
  H->rhs  = malloc(sizeof(fftw_complex)*((H->n)+18));
  H->rk   = malloc(sizeof(fftw_complex)*((H->n)+18));
  G->r    = malloc(sizeof(double)*((G->n)+18));
  H->r    = malloc(sizeof(double)*((H->n)+18));
  for (int jj = 0; jj < 7; jj++) {
    G->k[jj] = malloc(sizeof(fftw_complex)*((G->n)+18));
    H->k[jj] = malloc(sizeof(fftw_complex)*((H->n)+18));
  }
  fhlog = fopen("run.log","a");
  for (int k = 0; k < (G->n)+18; k++) { 
    G->r[k] = G->a + (k-9)*(G->h);
  }
  for (int k = 9; k < (G->n) + 9; k++) {
    if (k % 2 == 1) {
      G->psi[k] = LFT->psi[5 + k/2];
    } else {

      m0 = l0*(LFT->psi[k/2 - 5] + LFT->psi[k/2 + 14]);  
      m1 = l1*(LFT->psi[k/2 - 4] + LFT->psi[k/2 + 13]);
      m2 = l2*(LFT->psi[k/2 - 3] + LFT->psi[k/2 + 12]);
      m3 = l3*(LFT->psi[k/2 - 2] + LFT->psi[k/2 + 11]);
      m4 = l4*(LFT->psi[k/2 - 1] + LFT->psi[k/2 + 10]);
      m5 = l5*(LFT->psi[k/2 + 0] + LFT->psi[k/2 +  9]);
      m6 = l6*(LFT->psi[k/2 + 1] + LFT->psi[k/2 +  8]);
      m7 = l7*(LFT->psi[k/2 + 2] + LFT->psi[k/2 +  7]);
      m8 = l8*(LFT->psi[k/2 + 3] + LFT->psi[k/2 +  6]);
      m9 = l9*(LFT->psi[k/2 + 4] + LFT->psi[k/2 +  5]);

      /*m6 = t0*(LFT->psi[k/2 + 1] + LFT->psi[k/2 +  8]);  8th order
      m7 = t1*(LFT->psi[k/2 + 2] + LFT->psi[k/2 +  7]);
      m8 = t2*(LFT->psi[k/2 + 3] + LFT->psi[k/2 +  6]);
      m9 = t3*(LFT->psi[k/2 + 4] + LFT->psi[k/2 +  5]);*/

      m1 = m1 + m0; m2 = m2 + m1; m3 = m3 + m2; m4 = m4 + m3; m5 = m5 + m4; 
      m6 = m6 + m5; m7 = m7 + m6; m8 = m8 + m7; 
      G->psi[k] = m8 + m9;

      /* m4 = l4*((LFT->psi[k/2+2]) + (LFT->psi[k/2+3]));
      m3 = l3*((LFT->psi[k/2+1]) + (LFT->psi[k/2+4]));
      m2 = l2*((LFT->psi[k/2+0]) + (LFT->psi[k/2+5]));
      m1 = l1*((LFT->psi[k/2-1]) + (LFT->psi[k/2+6]));
      m0 = l0*((LFT->psi[k/2-2]) + (LFT->psi[k/2+7]));
      G->psi[k] = m0 + m1 + m2 + m3 + m4; */
    }
    fprintf(fhlog,"%.14e\t%.14e\t%.14e\n", G->r[k], creal(G->psi[k]),
                                  creal((p->a)*exp(-pow((G->r[k])/(p->sigma), 2)))  );
  }

  fprintf(fhlog,"\n");
  for (int k = 0; k < (H->n)+18; k++) {
    H->r[k] = (H->a) + (k-9)*(H->h);
    H->psi[k] = LFT->psi[(LFT->n)/2 + k];
  }
  dt = 0.25*dt; 
  free(LFT->psi); free(LFT->wrk); free(LFT->wrk2); 
  free(LFT->rhs); free(LFT->rk);  free(LFT->r);
  for (int jj = 0; jj < 7; jj++) free(LFT->k[jj]);
  free(LFT);
  LFT = G; nT = nT + (H->n)+18;
  p->left = G;
  update_boundaries();


  /*mgl_delete_data(xx);
  mgl_delete_data(yy);
  mgl_delete_data(zz);
  mgl_delete_data(aa);
  xx = mgl_create_data();
  yy = mgl_create_data();
  zz = mgl_create_data();
  aa = mgl_create_data();
  mgl_data_create(xx,nT,1,1);
  mgl_data_create(yy,nT,1,1);
  mgl_data_create(zz,nT,1,1);
  mgl_data_create(aa,nT,1,1);*/
}






